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Abstract 

A recurring task in image processing, approximation theory, and the 
numerical solution of partial differential equations is to reconstruct a 
piecewise-smooth real- valued function f{x), where x G R'^, from its trun- 
cated Fourier transform (its truncated spectrum). An essential step is 
edge detection for which a variety of one-dimensional schemes have been 
developed over the last few decades. Most higher-dimensional edge de- 
tection algorithms consist of applying one-dimensional detectors in each 
component direction in order to recover the locations in R'^ where f{x) 
is singular (the singular support). 

In this paper, we present a multidimensional algorithm which identifies 
the wavefront of a function from spectral data. The wavefront of / is the 
set of points {x,k) € R^ x {S^~^ /{±1}) which encode both the location 
of the singular points of a function and the orientation of the singularities. 
(Here 5*^"^ denotes the unit sphere in A'^ dimensions.) More precisely, k is 
the direction of the normal line to the curve or surface of discontinuity at x. 
Note that the singular support is simply the projection of the wavefront 
onto its x-component. In one dimension, the wavefront is a subset of 
R^ X (S^/jil}) — R, and it coincides with the singular support. In 
higher dimensions, geometry comes into play and they are distinct. We 
discuss the advantages of wavefront reconstruction and indicate how it 
can be used for segmentation in magnetic resonance imaging (MRI). 

1 Introduction 



Consider an image, i.e. a function p : ^ R+. If the image is smooth (e.g. 
C"*(IR.^)), then the Fourier transform of p{x), denoted p(fc), will decay rapidly 
(and hence be localized near fc = 0). Discontinuities in the image cause p{k) 
to decay more slowly as |fc| — s- oo. Thus, information about the discontinuities 
can be said to be encoded in the high frequency components of p{k). The goal 
of spectral edge detection is to recover the location of the discontinuities from 
limited (and often noisy) information about p{k). 

As an example, consider the function p{x) — 1b{x) which is equal to 1 inside 
B — {x : \x\ < 1} and elsewhere. The set of discontinuities of this function (the 
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singular support) is given by {x : \x\ — 1}. One natural approach to computing 
the curve on which the discontinuities lie is to first find a set of point which 
lie in the singular support, followed by an algorithm aimed at connecting these 
points sets into a finite number of curves. In our example, the output would be 
the unit circle. A relatively recent and important class of methods for locating 
the singular support is based on concentration kernels ^U\, [TH [13 [Hj , a high 
pass filtering approach, which we describe briefly below. 

For the function p{x) — Isix), the normal at each point of discontinuity 
X is simply the normal to the unit circle at that point. The wavefront of this 
function is, therefore, {(a;. A:)}, with |a;| — 1 and k \\ x. In this paper, we study 
the problem of extracting the wavefront from continuous spectral data available 
in a finite frequency range |fc| < kmax in two dimensions. This extra information 
is useful practically as well as theoretically. First, it is easier to reconstruct 
curves of discontinuity from points in the wavefront than points in the singular 
support, both in closing "gaps" and in associating points on close-to-touching 
curves to the correct one [H [3 [SJ (HI HI] . Second, the directional information is 
useful for noise removal. If spurious points are included in the wavefront set, the 
normal (or tangent) data allows us to filter it out; it is unlikely that a random 
point and a random tangent will be consistent with the points and tangents 
that come from the actual curves of discontinuity, as we have shown previously 
in [n]. 

Our approach to edge detection is based on applying concentration kernels 
(high pass filters) to angular slices of the Fourier data. Rather than recovering 
the points on the edges (as in [101 US EDj we also determine the direction 
of the normal. In section [21 we present a precise statement of the problem to 
be solved and an overview of the full edge detection procedure f section [2^ . 
In section [31 we present a detailed analysis of the asymptotic behavior of the 
Fourier transform of the characteristic function of a smooth region. Section [4] is 
devoted to a discussion of the directional filters used to extract wavefront data 
and Section [S] describes the full algorithm. Many of the proofs are technical and 
we have relegated most proofs to the appendices. We discuss the application of 
our method to magnetic resonance imaging (MRI), where raw data is acquired 
in the Fourier domain, extending the algorithm of [H]). Finally, we note that 
our algorithm is closely related to the recent paper [18] , which describes a wave- 
front extraction procedure based on the shearlet transform, extending the use 
of curvelets in [5l [6] . 

2 Mathematical Formulation 

Let p{x) be a piecewise smooth function form M, and let xi^) be a C°° 

compactly supported radial function. The singular support of p{x) consists of 
the points xq G for which xiM^ ~ xo))f{x) has slowly decaying Fourier 
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transform for every A > 0, i.e.: 



VA > 0, sup 

\k\>kr 



e p{x)x{{x - xo)X)dx 



(2.1) 



In order to find the singular support, concentration kernel methods [10|, 
[El im im US] multiply the Fourier data p{k) by a function which gives heavier 
weight to high frequencies than to low frequencies (a high-pass filter). Since high 
frequencies encode the location of singularities but are unaffected by smooth 
parts of the image, this method isolates discontinuities from the rest of the im- 
age. In short, concentration kernel methods find the location of singularities by 
flagging local maxima in the inverse Fourier transform of the high-pass filtered 
Fourier data. 



The wavefront of a function consists of the points (xo,fco) G 



X Si for 



which the Fourier transform of xi'^i^ ~ ^o))pix) decays slowly in the direction 
ko = {kr, ke): 



VA > 0, sup 



.jrko-^^ p{^x)x{{x - Xo)X)dx 



= 0{K 



-3/2N 



(2.2) 



As indicated in the introduction, while the singular support of p{x) only contains 
the location of singularities, the wavefront also contains the direction of the 
singularities. 

Remark 2.1 In the language of computational geometry, the singular support 
is a set of points, while the wavefront is a set of surfels (pairs of the form (a;, k) 
with X representing a position and k a direction). 



2.1 Definition of the image class 



To simplify the theory, we consider a special class of images. In particular, we 
consider two-dimensional images supported on [0,1]^ and vanishing near the 
boundaries, which consist of a set of piecewise constant functions on which is 
superimposed a globally smooth function: 



p{x) 



M-l 



(2.3) 



where jj(t) are simple closed curves, and 17^(2;) = 1 for x in the interior of 7j 
and elsewhere. The "texture" term ptcx(a^) is band limited, i.e. Ptex(fc) — 
for k > fctex- 

Definition 2.2 Let^j{t) he a simple closed curve. The curvature at each point 
is denoted by Kj{t), the normal to ^j{t) is denoted by Nj{t), and the tangent is 
denoted by Tj (t) . 
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Figure 1: An illustration of Assumption [2l The black arrow illustrates the 
condition (|2.5ap , while the red arrow illustrates the condition (j2.5bp . 



Assumption 1 We assume the curves have bounded curvature, i.e.: 

Assumption 2 We also assume that the curves are separated from each other 
(see Fig. Qp.- 

svip\^,{t)--fj{t')\>S fori^j (2.5a) 
t,t' 

and that different areas of the same curve are separated from each other, i.e. if 
the curves are parameterized to move at unit speed, then: 

sup hdt)-j,{t')\>S (2.5b) 

|t-t'|>S-i7r/2 

In order to extract the wavefront, we assume each edge is associated with 
a sufficiently large discontinuity. Since our wavefront detector decays some- 
what slowly away from the wavefront, we also assume that the discontinuity is 
bounded from above, so as not to pollute nearby edges of lower contrast. We 
formalize this as: 

Assumption 3 We assume the contrast of the discontinuities in the image is 
bounded above and below: 

< P < Pj <P (2.6) 
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The data we are given and wish to segment are noisy samples of p{k) ob- 
tained on a 2™ X 2™ grid with spacing 27r in fc-space centered at the origin: 
k e 27rZ2™ with = {-2™"!, . . . , 2™"^ - l}^. Our segmentation goal is to 
recover the curves 7^ (t) . 

With the preceding sampling, we may define the maximum frequency avail- 
able in the image as /cmax = 2tt2"^~^. We assume that fctex < ^max- More pre- 
cisely, we assume that fcmax — ^tcx ^ 12 • 27r, providing at least 12 lattice points 
in the sampling beyond fctcx- In most of our experiments, we take /ctcx = 27r • 16, 
and fcmax — 27r • 32. 

Finally, we make the technical assumption that the curves have non- vanishing 
curvature: 

Assumption 4 We assume that the curvature of the curves is bounded below. 

Vi = . . . M - 1, K,{s) > K > (2.7) 

Remark 2.3 The assumption (j2.7|) implies that the region bounded by ji{t) is 
convex. This geometric fact is not used by our algorithm in any way. 

Remark 2.4 Assumption |4] is introduced only because it is required below 
by our proof technique for the correctness of the directional filters. It is not 
strictly necessary, and it would be straightforward to extend our analysis to cases 
where the curvature vanishes. This, however, would require more complicated 
conditions on higher derivatives ofjjit) in places where the curvature vanishes, 
and correspondingly more complicated proofs. See Figure [6] which illustrates 
that our algorithm works even when Assumption [4] is violated. 

2.2 Wavefront Extraction Methodology 

The algorithm we present in this work takes an image, given as spectral data 
p{k), and extracts a set of surfels sampled from the wavefront. We do this in 
two steps. 

First, we construct a set of directional filters: 

[De,c.p]ix) = T-'[W{k,)V{ke - e)p{k)] (2.8) 

Here, VV(fcr) is a high pass filter in the radial direction (the radial filter), V is 
a smooth function, compactly supported on [—a, a] (the angular filter), and T 
is the Fourier transform. Note that, given a direction 6, the angular filter is 
supported on the angular window [9 — a, 9 -\- a]. The angular and radial filters 
will be related by the parabolic scaling: 

width'^ = length 

(Equivalently, they will be chosen to satisfy width = length^ in the ai-domain.) 
In particular, this implies that the filter angle a ~ (pass-band) where 
pass-band is the center of the pass-band of the radial filter (to be determined 
in section [4|). 
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When applied to the function p{k) , the directional filters will return a func- 
tion which is small except near locations where Nj{t) points in the direction 9. 
This allows us to pinpoint the locations of singularities with direction 9. The re- 
sult is, in some sense, a directional version of the jump function of pUllT^fT^fTB] . 
Spikes (local maxima) obtained from these directional filters correspond to sur- 
fels in the wavefront of the image. 

Remark 2.5 For the algorithm described here to work, it is required that ac- 
curate values of the continuous spectral data be available. The discrete Fourier 
transform (DFT) of p{x) can not be used as a substitute for the continuous 
data, since aliasing induced by the DFT will destroy the asymptotic expansion 
of p{k). 

2.3 A note on our phantom 

In the medical imaging literature, the Shepp-Logan phantom (which is piecewise 
constant) is a traditional choice for analysis and validation purposes. We have 
added a smoothly varying component (see (|2.3p ). with ptcx(a;) defined as a sum 
of Gaussians whose bandwidth (to six digits of accuracy) is equal to half the 
bandwidth of the measurement data. 

2.4 Informal Description of the Algorithm 

The full algorithm proceeds in three steps. 

1. Construct directional filters, based on fcmax (the maximum frequency con- 
tent of the data) and fctcx (the frequency content of the "smooth" part of 
the image). 

2. Multiply the data p by each directional filter indexed hy 9 — jir/A, j = 
1, . . . , A) and transform back to the image domain. 

• Apply a threshold in the image domain. 

• Add each point x above the threshold to the surfel set as the pair 
{x, ke). 

3. Given the set of all surfels, use the algorithm of ,T7] to reconstruct the set 
of curves 7j(t) that define the discontinuities. 

The last step is described briefly in section [6l while the bulk of the paper is 
devoted to the surfel extraction procedure itself. 

3 Large k Asymptotics of l-yjik) 

In this section we wish to compute the large k asymptotics of 1^^. (fc). In partic- 
ular, we will show that ijj{k) is dominated solely by the parts of 7j(i) where 
Y^{t) • fc = 0, i.e. where k _L Tj{t). 
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Figure 2: A plot of the phantom used in this work. Tlie image in the sec- 
ond panel is obtained by applying the DFT to the (exactly known) continuous 
Fourier Transform of the phantom on a 64 x 64 grid of samples in fc-space. 
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We can use Green's theorem to rewrite the Fourier transform of l.y. (x) 

-2 



as follows. Let F{k,x) = {Fi{k, x), F2{k, x)) 
[—k2,ki]'^. Then by Green's theorem: 



s*'='"=fc-'- with 



l7.(fc) 



dxiF2{k,x) - dx2Fi{k,x)dxidx2 



FikMt))-'^dt^ ' 



dt 



'fe-7.W^i . V.(t)dt (3.1) 



(with the region bounded by 7^). This trick is taken from [22l. 

We can now use stationary phase to approximate ly.{k) for large k. For 
this, express k in polar coordinates {kr,kg), fix a direction kg, and consider 
what happens as kr becomes large. As we remarked earlier, the phase k ■ ^j{t) 
becomes stationary only when k ■ J jit) = or fc • Tj{t) = 0. This is precisely 
where k points normal to the curve, and it is these locations that dominate 
l,^(fc): 

Proposition 3.1 Let tj{k) correspond to the value oft at which Nj{t) || k and 
Nj{t) ■ k > (i.e. the normal to 7j(t) points in the direction k). Then: 



Af-l Af-l 



\k 



3/2 



where E{k) < <Lgeo with 



Mp 4 



Kj{tj{k)) 



3/2 



3p 



sup, II 7'" (i) I 



i^jitji-k)) 
E{k) 



(3.2) 



■ ^ arclength(7j ) . (3.3) 



Note that £goo incorporates both geometric and contrast information about 
the image itself. This result is proved carefully in Appendix [Xj The basic idea 
behind the proof is simple, however. Set k = k^kg, with kg fixed. Then: 



e»^-7.(t)fci . ^'{t)dt = / e'^^-^C^fc-^ • Mt)dt 



with fit) — kg ■ J jit). The phase function fit) is stationary when kg ■ J jit) = 0, 
or equivalently the place where kg ■ Njit) — ±1 (i.e. tj(fc) = tjikg)). By 
Assumption m we find that Jjit) — Kjit)Njit) is nonzero. We restrict the 
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consideration to an interval I = [tj(kg) — e, [tj(kg) + e] and apply stationary 
phase: 




,tKf{t,(ko))^i^ . j'^{tj{kg)) + remainder 



= / !l^^e*^-''^"(*^"(^» + remainder 

Proving Proposition [3Tl] is done by adding this result up over all the curves, and 
all points of stationary phase, and estimating the remainder. 



3.1 WhatifK = 0? 

It is important to note that all is not lost when k = 0. In this case, although 
the coefficient on k^'^/'^ in (|3.2p becomes singular, the asymptotics of 1-y. (k) do 
not. 

In this case, what happens is that the leading order behavior becomes 
0(fc^"^"^/"), where n is the order of the first non- vanishing derivative. This can 
be seen relatively easily from stationary phase, although rigorous justification 
requires a long calculation. However, the limiting case n = cxi (a straight line) 
is easy enough to treat. 

Proposition 3.2 Let -f(t) ^ a + (b - d)t for t G [0, 1]. Then: 



e^-^^i)k^ . ^' {t)dt 



-1 k^-{b-a) 
jlp' fc-(b-a) 



k-{b-a)=0 



0{kr~'), k-{b-a)^0 
e^^-^0{kr~^), k-{b-a) = ^ ' ' 



Proof. If fc • (6 - a) 7^ 0, then: 

e'^^-^W^J- . y(i)rft = f e'''<(^-^')^+^'>k^ ■ {b-a)dt 

^ik-{h — a) 







k^-{b-d)e'^-^ / e'''-^^''~^'^'Ut = k^-{b-a)e'^-^ 
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ik ■ (b — a) ik ■ [b — a) 
k^ ■ {b-d) 



ik ■ (b — d) 
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Multiplying hy (i k ) ^ yields the result we seek. 

The asymptotics are straightforward to compute (from the second line of 
(|3.4p ). but note that the constant in the 0{kr~'^) term in (|3.4|) is not uniform 
in kZ{b- a). □ 



4 Directional Filters 

We are now in a position to build the filter operators Dg^a of (|2.8p . which will 
allow us to extract edge information from the signal. We demand that the radial 
filter takes the form 

Wikr) = W^ikr - (fcmax + fctox)/2) + (kr + (fc„,ax + A:tcx)/2) (4.1a) 

where yVP(r) is a positive, symmetric function. This means that the 1-dimensional 
inverse Fourier transform of yV(fcr) is 

>V(r) = cos j^ ^'°^- + ^tex ^^ ^p^^^ ^^^^^^ 

The intuition behind this operator is the following. Multiplying by yV{kr)V{kg — 
0) localizes on the region k « {kr cos 6, kr sin 9). Dropping all but one of the 
terms in Proposition [3Tll we find that: 

yV{kr)V{ke - e)Q{k) « e'''-^^^'^^^^\j{tj{k))-^'^A[k) . 

The term A{k) incorporates both the l/|fcp/^ decay from the asymptotics of 
l^j. {k) and the localization to ke and large kr from the filters. 
Therefore, if we inverse Fourier transform, we will obtain 

Provided A{x) is a sharply localized bump function, this will be a bump located 
at ^j{tj{k)). Of course, this calculation is not exactly correct, and is presented 
merely to obtain intuition. We will go through the details shortly, but require 
a few definitions. 

Definition 4.1 Given W^k,.), define the auxiliary functions: 

W(r) = J e-''^^-W{kr)dkr (4.2a) 
W(i?) = sup|>V(r)| (4.2b) 

r>R 

Definition 4.2 Define the set Sle.a C to he the set of points where some 
Jj{t) has normal pointing in the direction [9 — a, 9 + a], i.e.: 

M-l 

2le,a = IJ [{-fj{t)}te[tjie-a).t,ie+a)]^ {lj{t)}tGitj{9-a+7T),tj{e+a+7i)]] (4.3a) 

3=0 
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We also define ^g^a to he the set of arcs excluding ^j(tj{9)). 



^e,a = [J [{ji{t)}telUie-a),U{e+a)] U {'-ft{t)}t^[ti{e-a+-R}MS+a+7^)]] (4-3b) 

The goal of our directional filters is to approximate the location of Slg.Q. That 
is to say we want Dg_aP to be large near Slg^Q and small away from it. The 
decay in the tangential direction away from a point in Sle^a is at least of the 
order 0{r-^). To ensure that the decay in the normal direction is as fast, we 
consider the case when >V(r) — 0{r^^) (see Fig. [3] d). 

We then have the following result which proves the directional filters ap- 
proximate the location of Slg.Q. 

Theorem 4.3 Suppose that a satisfies 



(4.4a) 



and W{kr) satisfies (|4.1ap as well as the following decay condition: 

W(r) < ^ 
r 

Then away from Slg^Q, the directional filter has the following decay: 

G:(W,V,a) 



(4.4b) 



M-1 
3=0 



< 



d{x,^g^a) 



\Vikg)\\ 



(4.4c) 



where 



2MpV2TT 



2 ||VV(fcr)/fcr||^i(][j ) 



llv'(MllLi + llv(fc«)IL. 



|7'"W| 



2^2 



+ ^\\W{z){z+l)\ 



(4.4d) 



At the point x = Jj{tj{0)) 

[9 ■ Nj{t)]Dg^^ 



AI-l 
3=0 



> T(W, V,a) 



(4.4e) 
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with 



1(}V,V,a) = J —p inf Wf(r) 

V 2K-re[-a2/2K,Q2/2K] 



g:(W, V,a)(2M- 1) 
2MS 



-1/2 



(4.4f) 



We postpone the proof of this result to Appendix [Bl where we compute the 
leading order asymptotics of a single segment of the curve and put the various 
pieces together. 

Remark 4.4 While the expressions above are somewhat involved, for the filters 
described in the next section, (|4.4cp simplifies to 



M-l 



< 



d(a;,2te,a) 



O 



1 



1/2, 



1/2 



(4.5) 



Note that the second term in this estimate becomes negligible as A:niax increases. 
The term £(W, V, a) determines the rate of decay of the directional filter away 
from a surfel and it should be as small as possible (to minimize noise), for which 
the norm of W should be small. At the same time, we want ICW, V, a) to be 
as large as possible (to maximize signal). For this, the norm of W should be 
large. We will need to balance this competition. 

Remark 4.5 If the decay of >V(r) is faster than (j4.4bp . we can of course get 
sharper results. However, this would require extra geometric conditions and 
would require Theorem 14.31 to make distinctions concerning the direction from 
X to Stg.Q. To simplify the exposition of this paper, such considerations will be 
reported at a later date. 



4.1 The estimate is suboptimal 

Let us consider the application of our directional filters on a simple numerical 
example. We consider a 64 x 64 pixel image on [0,1]^. The image is taken 
to be the phantom described in Section 12.31 and the parameters are /cmax = 
647r « 201.1, /fctex = 327r « 100.5, a ^ tt/16, k = 0.1 and ||7"'(i)||^^ < 5. 
With this set of parameters, £(yV, V, a) — 0.3 if we use the windows V{k0) = 
{2a)~''-l[Q_a.e+a]ike) an d W(fc r) = (fcmax-fctcx)""^l[fet<3^,fc„^^](fcr). The first term 
on the right-hand side in (I4.4cp is, therefore, approximately 0.3/(1/64) « 25 if we 
move one pixel away from a surfel. The second term, however, is approximately 
85 so that we have no reason to expect that the directional filters will yield any 
useful information. However, numerical experiments show that they do in fact 
work even in this case. Figure [3] shows that the directional filters do yield the 
location of the edges. Figure H] shows that even in the presence of noise (up to 
7.5% of the total image energy), the directional filters still yield correct results. 

It is also useful to compare the directional filters to a naive algorithm, namely 
computing the directional derivative. As is apparent from Figure [5l this naive 
method does not perform as well as our algorithm. 
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Figure 3: An illustration of the sharp directional filter. The first panel shows 
the image. The second panel shows the directional filter in fc-space. The third 
panel shows the edge map [Dg^aP]{x) with 9 = n/A and a = 7r/16. The red 
lines are the | [Dg^aP] — 2.4 contour lines, while the green lines are the actual 
(analytically known) edges of the image. The fourth panel shows the directional 
filter in the x-domain. 
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Edgemap, noise = 2.5% 
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Edgemap, noise = 7.5% 
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Figure 4: An illustration of the directional filters in the presence of noise. The 
image and parameters are the same as in Figure[3l The figures in the left column 
are plots of [Dg^aP]{x)- The gray dotted lines in the right column illustrate the 
discontinuities of the image (which we know exactly, c.f. Section l2.3p . while 
the red regions are the places where |[D6i_ctp](a;)| = 2.4. Noise amplitude is 
described as a percentage of the total image strength. Noise begins to corrupt 
the data when it reaches 7.5%, and becomes a serious problem at 10%. 
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Directional Derivative 
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Directional Derivative, Thresholded 
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Figure 5: Result of edge detection based on measuring the directional derivative, 
i.e. djr/4:P{x). The green lines represent the analytically known edges of the 
image, while the red regions are the regions where |97r/4P(a;)| > 150. 
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Exact Contours, Extracted Surfels 
1.0| i i i i — 




Q^- IIMtltlllltltl 



0.2 



O.Ql ' ' ' ' 1 

0.0 0.2 0.4 0.6 0.8 1.0 



Figure 6: Result of edge detection applied to an image where Assumption |4] is 
violated. The green lines in the second panel indicate the actual discontinuities, 
while the red arrows indicate extracted surfels (at angles 9 = and 9 = tt/A. 
The vertical and horizontal Gray bars are artifacts from using the DFT to 
reconstruct a discontinuous image. 
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4.2 Parabolic Scaling 

The choice of a is an important one. We want a to be as small as possible, 
since this will give us better angular resolution. 

On the other hand, the constant bounding the size of the filtered image 
away from the edges is directly proportional to oT^ ^ as we will show shortly. 
For simplicity, let us impose the constraint that 



l^(fce)llLi(fe<,,dfe, 



1. 



(4.6) 



Now let us take V(fce) to be a fixed filter scaled with a, i.e. V(/ce) = oT^VikQ /a), 
where V[kg) is supported in [—1,1] and |l^(^e)|lLi(d*:e) ^ ^- With this choice 
of V(fce), we find that: 

l|V'(^e)||Li(si,dfe,) = ||a-^y'(fce/a)||^i(gi ^^^^ = ll^'(fce/a)llLi(si,dfee) 

(4.7) 

Substituting this into (I4.4dp and assuming a to be very small yields: 



e:(w,v,a) 



'2ti 



r'(fc./a)IL.(s.,,,,) + |lV(fc,)|l^, 



2k2 



+ ||W(^)(^ + ) = O («;-i/2^-i ||W(fc.)/fc.|l^,(^ ,,^)) (4.8) 



To prevent £(W, V, a) from being too large, we want to ensure that a is as large 
as possible. 

However, to ensure that the filter is sufficiently large on the edges, we need 
to prevent a from being too large (Theorem 14. 3p . In particular, we require (c.f. 

6311) 



a < 



(4.9) 



To satisfy both these constraints, we simply replace the "<" in (|4.9p by ' 
This yields the standard parabolic scaling used elsewhere O [6l [HI [211 [22] . 



Remark 4.6 The requirement ()4.4a|) yields the standard parabolic scaling used 
to analyze line discontinuities in harmonic analysis O [6l [Til [211 123] . In the x- 

domain, the standard parabolic scaling uses elements with width ^ length^. In 
the /c-domain, this translates to width ^ ^length. 



This also implies that for the filter to approach zero as fctox, k„ 



oo, we 



require that || W(fcr)/fcr|lz,i( 



o(kt. 



-1/2 
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4.3 Choosing the Filters 

We now consider the simplest choices of window possible which satisfy our 
assumptions. We take W(fcr) to be a a step function, i.e.: 

1 



W{kr) = 



(l[/ctex,fc„,,x](fcr) + l[-fe.„,x,-fct„,](fcr)) 



2(^max ^tcx 

and we take V{ke) to be a triangle function: 

V{ke) = -\max{a~ \kg \ , 0} 
In this case, it is straightforward to show: 

kr'^/^Wikr) 



|ww(- + i)ILo.(«,,.) 

ll"^'(fce)|lii(si,dfe«) 
l|V'(fc«)ILHs^dfe«) 



"'max r A,tcx 
ln(fcmax/ fctcx) 



< 



^max 


" ^tcx 


1 




h — 
"-max 


^tcx 


1 






^tcx 


1 





Then we have: 



^COs(2a)(fcmax " fctox) 
2 ln(fcmax/fctox) 2 



= 1 



= 2a-^ = 2 



max < 1 , 



max T' n^tcx 



fcto 



2«;2 



Similarly, the remainder term has the bound: 

1 



n ^^^^ 



< 



1. 1/2 ,L 1/2 



|v(fc,)|| 



(4.10) 



X £, 



gco , 



where Cgeo is defined in p.3p 
Therefore, we find that: 





M-l 

















o 



ln(fcmax/fctox)\/fc max ~r ft-tcx 



fc. 



max "-tcx 



U 1/2 , K 1/2 
"-max ~ "-tcx ' 



(4.11) 
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Thus, away from Stg^Q, the directional filter approaches zero. The 0(ln(fcmax/^tcx)/(fci 
fctcx)) factor is present because we are attempting to extract spatial information 
from a frequency band of width /cmax — ^tox- The \/ {kms,^^'^ + ^tex"'^''^) term 
is present because the asymptotic expansion we used to derive the directional 

I/O 

filters decays only kr ' faster than the leading order terms. 



5 Surfel Extraction 

Theorem 14.31 proves that, provided we choose the parameters correctly, direc- 
tional filters will decay away from 2le^Q,. Extracting surfcls from the filters is 
therefore simply a matter of choosing the parameters correctly and seeking local 
maxima. 

We know that at the point '^j{tj{0)), (|4.4ep provides a lower bound on the 
size of the directionally filtered image. We also know that away from Slg^cn 
(|4.4c[) provides an upper bound on the size of the filtered image. 

We wish to say that if the filtered image is "large", then we are near Slg.Q, 
otherwise we are not. Therefore, we will need the lower bound at ^j{tj{0)) to 
be greater than the upper bound away from SI^.q: 



TT 



2R 



Pmke)\\LH,u,) ^, inf ^^^W^r) 



g:(>V,V,a)(2A/- 1) 
2Md{x,'QieJ) 



L^{R,dkr.) 



\Vike)\\ 



^gco 



> 



£(>V,V,q) 



+ 



kr-^/'^W{kr) 



Li(R,dS,v) 



Cgeo (5.1) 



Let us also use the convention that 



l|V(fc«)ILi(d..)-l- (5.2) 
Using the fact that d{x,'Qig^a'') > d{x,^g^a), we find that (|5.1|) implies: 



pj— inf wp(r)/e:„e 

— V 2k rel-a^/2K,a^/2K] 



^ (4M-l)g:(>V,V,a) ^ 
2Md{x,^e,a) 



-1/2 



Wikr) 



L^(R,dkr.) 



(5.3) 



Thus, if kr ^^^W(kr) is sufficiently small, then (15. ip will be true 

Li(K,dfer) 

whenever d(x,^g_a) is sufficiently large. In particular, if we make each term 
on the right side of (|5.3p smaller than half the left side, this equation will be 
satisfied. 
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We summarize this result in the following corollary to Theorem I4.3[ which 
shows that the directional filter is large only when d{x, ^e,a) is sufhciently small. 



Corollary 5.1 Suppose that: 



Li(R,dfc,.) 

Define the surfel location error to be 



< \pJ^ inf yV^{r)/<Lneo. (5.4) 

4-V 2k re[-aV2«;,aV2Kl 



S(>V,V,a) 

_ 2 f2R 

p\ IT 

Then whenever 
we have that 



mi WP(r) 

rG[~a^ / 2k, /2k] 



{AM - 1)€{W, V, a) 



2M 



-geo ■ 



d(x,2te,„) >D(>V,V,a), 



De.a 



Af-l 
3=0 



(x) 



< s:(w,v,a). 



(5.5a) 
(5.6) 
(5.7) 



Thus, when x is located a distance at least S(yV,V,Q;) from Sle.a, the filtered 
image is smaller than 1{W,V,a). On the other hand, at the point jj{tj{a)), 
we know that the filtered image is larger than T(yV, V, a). Thus, we obtain the 
following thresholding algorithm for locating surfels in the wavefront: 

Algorithm 1: Surfel Extraction 

Input: The image in the Fourier domain, i.e. p{k) and a desired minimal 
sampling rate e. 

Output: Surfels which approximate the wavefront of p{x) in the 
direction 6. 

1 let fi{x) := [Dg^aP]{x). 

2 let Z:={xe : |/(x)| > T(>V,V,a)}. 

3 Cluster the set Z. Any two points are part of the same cluster if they are 
located a distance D{W, V, a) apart. Let S denote the set of clusters. 

4 let RESULT :== [] (the empty set) 

5 foreach s £ S do 

6 Let the midline of s be the set {midpoint((a; + 9M.) D s) : x E s}. 

7 Sample the midline of s with spacing at least e, calling the result Q. 

8 foreach q E Q do 

9 I Add the surfel (q, 9) to RESULT. 

10 end 

11 end 

12 return RESULT 
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Figure 7: An illustration of the result of Algorithm [T] The red arrows indicate 
the location and direction of extracted surfels in the direction 7r/4 (with a = 
7r/16). The surfels are overlayed on the image (reconstructed by DFT). 

An example of Algorithm [1] applied to the same image as in Section 14.11 is 
shown in Figure [T] The algorithm generates no surfel a distance more than 
1.0/64 (i.e., one pixel) away from the actual edge. 

We have the following result concerning correctness of Algorithm [T] 

Theorem 5.2 Suppose that in addition to (|4.1ap . (|4.4ap . (|4.4bp . (|0)) . 

the following constraint is satisfied: 

S(W,V,a)<- (5.8) 

Then for every {x,6) in the result of Algorithm]^ there is a corresponding surfel 
{x',0') in the wavefront of p{x) with the property that: 

\x- x'\<D{W,V,a) (5.9a) 

\e~e'\<a (5.9b) 

Additionally, for each j , at least one surfel output by Algorithm [7] will ap- 
proximate some surfel {jj{tj{X)),Nj{tj{X))) in the arc segment arc segment 
{j,it,{X)):\e[e-a,e + a]}. 
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Proof. Define the arc segment: 



= {ij{t.j{\)) : A e [61 - a, e + a]} 

By Corollarv l5.ll any point x £ Z is located a distance at most !D(yV, V, a) 
away from %e,a- 

Now, consider any segment s oi Z . For any a; e s, there is some point x' in 
some arc segment Aj for which d{x,x') < !D(>V, V,Q!). We need to show that 
all points y € s are located at most a distance J)(W,V,a) from the same arc 
Aj. Consider a point y with d{y,Ak) < ©(W, V,q;). Suppose also that y' e A^ 
is a point for which d(y, y') < Sur f el Resolution. Then: 

S < d{x',y') < d{x',x)+d{x,y)+d{y,y') < D{W,V,a)+d{x,y)+D{W,V,a) 

Subtracting 2!D(yV, V, a) from both sides and applying (|5.8p imphes that d{x, y) > 
J/3 > J)(W,V,Q!) and therefore x and y are not in the same segment. Thus, s 
consists only of points a distance S)(>V, V, a) from Aj. 

Therefore, any point on the midline of S is located a distance at most 
V, a) from SI^.q. By definition, any point on Sle^ct has a normal pointing 
in some direction in [6' — a, + This proves the Surfels returned by Algorithm 
[1] accurately approximate surfels in the wavefront of p{x). 

To prove that each arc segment has at least one surfel in it, note that by 
Theorem 14. 3[ the point ^j{tj{9)) is contained in Z (c.f. (|4.4cp V This implies 
that each arc segment {7j(tj(A)) : A G [6* — a, 6* + a]} generates at least one 
segment s of Z. There will be at least one sample taken from this segment, 
which will generate a surfel in the output of Algorithm [1] Thus, Theorem 15.21 
is proved. □ 



6 Segmentation: connecting the surfels is better 
than connecting the dots 

As we have indicated earlier, one of the reasons for developing a surfel/ wavefront 
extraction procedure is segmentation - by which we mean the reconstruction of 
the curves of discontinuity 7^ which divide the image into well-defined geometric 
sub-regions. 

The first step in reconstructing the curves is to reconstruct their topology. 

Definition 6.1 A polygonalization of a figure {'yj(t),j = 0, . . . , Af — 1} is a 
planar graph T = (V, E) with the property that each vertex p € V is a point on 
some "fj{t), and each edge connects points which are adjacent samples of some 
curve Jj{t) (see Fig. 0). 

There is a substantial literature in computational geometry discussing the 
task of taking as input a set of unordered points that lie on a set of curves and 
returning a polygonalization [51 [71 [HI IS] • In [H] ; we described an algorithm for 
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Curves and polygonalizations 



- - 7 




Figure 8: A curve and it's polygonalization. 

polygonalization that uses both point and tangent data (i.e. surfels) and showed 
that the method is significantly more robust. It is easier to remove spurious data 
from the set of surfels and the sampling requirements are much weaker (see Fig. 

m ■ 

In particular, Theorem 2.8 of [T7] shows that, given a set of points and tan- 
gents (a discrete sampling of the wavefront of p{x)), then there is an algorithm 
that returns the correct polygonalization provided S > 2Ke^ and e < (\/2k)~^, 
where e is the maximal distance between samples on the curve. 

If the point and tangent data are corrupted by noise (as they are in practice), 
then we need to assume a maximal sampling rate as well. Otherwise noise could 
change the order of samples on the curve. The following theorem provides 
technical conditions under which one can prove that the algorithm is correct 
when applied to noisy data. 

Theorem 6.2 f_?7 | Theorem 3.2] Suppose that Assumptions[l\ and[^ hold, that 
noise in the point data is bounded by and that noise in the tangent data is 
bounded by ^. Suppose further that: 

S > AC + AeC + 2.lKe^ , (6.1a) 

e<T^, (6.1b) 
KV2 

and that adjacent points on a curve are separated by a distance greater than 
[(1 + 2'^/^)(2^e + C)]. Then there is an algorithm that correctly reconstructs the 
figure. 
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Figure 9: The advantages of connecting surfels over connecting points are il- 
lustrated in this figure. First, the separation between points (e) can be much 
greater than the separation between curves (S) . Without information about the 
tangent, e and 6 must be of the same order. Using the algorithm of [T7], it is 
easy to automatically assign points to the correct curve. Note also that points 
such as C are easy to filter away when tangent information is available, even in 
the presence of modest amounts of noise. 



Once the polygonalization of the curve set has been obtained, one can ap- 
proximate the geometry with higher order accuracy. This is particularly easy 
in the case of surfel data; between each pair of points, cubic Hermite interpola- 
tion constructs a fourth order polynomial in arclength that interpolates the two 
points and matches the derivative (tangent) data as well. This achieves fourth 
order accuracy. 

The full segmentation algorithm follows. 

Algorithm 2: Segmentation 
Input: The Fourier transform of the image, p{k). 
Output: A set of curves approximating the discontinuities of p{x). 

1 let S = [] 

2 for j ^ 1 . . . A do 

3 let 6* jn/A 

4 let s — result of applying Algorithm [T] (the Surfel Extraction 
algorithm, see p. [20]) to p{k) in the direction 6. 

5 Append s to S. 

6 end 

/* Now S contains surfels pointing in the direction 6 for many values of 

e. V 

7 Call the algorithm of [17] to polygonalize S. 

8 return C, where C is the Hermite interpolant of the polygonalization of 
S. 



This algorithm can be proven "correct" in the sense that, for sufficiently 
large k^a,^, the algorithm will return a set of curves which are topologically 
correct. To do this, we need need to prove that the output of Algorithm [1] 
meets the requirements of Theorem 16.21 This requires verifying that (|6.1[) are 
satisfied. Note that we take Assumptions [l] [2l [3] and [4] as given. 

After extracting surfels from the image (as per the loop in lines 2-6 of Algo- 
rithm [2|), we find that the error in each surfel's position is bounded by Theorem 



24 




max "'tcx 

Th — ^3/2 
V^max "'tcx J 

1 

V^max + k 



tcx 



(6.2) 
(6.3) 



Both these quantities are 0(fcinax~^^^), holding all other factors fixed, and can 
therefore be made as small as desired. 

Note that the angle between adjacent surfcls returned by separate appli- 
cations of Algorithm [1] is at most 2a + tt/A. By Lemma [A. 21 the separation 
between two such surfels in arc length is at most (2a + tt/A)/k. By taking 
A = O(Vfcmax) (e.g. A = TT/a), we find that e = 0(fcmax"^^^) and thus (|6.1bp 
is satisfied for sufficiently large fc,nax- To satisfy (|6.1ap . observe that: 

4C + 4ee + 2.lKe2 < 0(fc„,ax"'/' + fcmax"' + A:„,ax"') - 0(fcmax"'/') < S (6.4) 

For sufficiently large fcmax, this implies (j6.f ap is satisfied. Thus, we have shown 
that (jO.fbp is satisfied. This implies that for sufficiently large fcmax (holding 
all other parameters fixed) , if Assumptions [T] and [2l [3] and [4] are satisfied, then 
Algorithm [2] will successfully segment the image. 
We summarize this in a Theorem. 
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Segmented Image 




Figure 10: The result of applying the segmentation algorithm to a 64 x 64 grid 
of spectral data. Note that the separation of two curves near the center is 
smaller than a pixel, but that wavefront/surfel reconstruction has no difficulty 
in resolving the them. 

Theorem 6.3 Suppose that AssumptionUl andj^ hold. Then for suffi- 
ciently large ktex, kmax, Algorithm\^ will successfully approximate the singular 
support of p(x) . 

The result of applying Algorithm [2] to spectral data for our phantom on 
a 64 X 64 grid is shown in Fig. [TOl The deviation from the exact result is 
noticeable, but is on the order of a pixel since we are using low resolution data 
to make the nature of the error clear. 
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7 Conclusions 



In this paper, we have described a new method for edge detection that can be 
viewed as an extension of the method of concentration kernels [TOl [T2l [T4l [T5l [T6] . 
We use more complex filters in order to recover information about the wavefront 
of a two-dimensional image rather than just its singular set. That is, instead of 
trying to locate a set of points that lie on curves of discontinuity, we look for 
both those point locations and the normal (or tangent) directions there. This 
allows us to reconstruct edges more faithfully and robustly, using the algorithm 
developed previously in [171 . We have focused here on a rigorous mathematical 
foundation for the method, based on detailed asymptotics and Fourier analy- 
sis. Although in this work we require that the curvature of the discontinuities 
not vanish, this assumption is merely technical. The algorithm works properly 
even when that assumption is violated, and even for singularities which are not 
differentiable (see Figure [H] for an example). A major advantage of extracting 
surfel information is that one can more easily "denoise" the data, as discussed 
in detail in l7] and illustrated in Fig. [5) A number of improvements can still 
be made, including the incorporation of nonlinear "limiters" to reduce the os- 
cillations produced in the physical domain from our linear filtering procedure 
(see, for example, J6,). 

Recovering local information about a function from partial Fourier data is a 
rather subtle issue, as demonstrated by Pinsky 20J who showed that spherical 
partial Fourier integrals do not converge pointwise to the characteristic func- 
tion of the unit ball in M.^. His analysis suggests that radial variations of the 
concentration method may not converge either (though of course appropriately 
filtered versions will). 

A limitation of the method described here is that we have assumed the 
image consists of a globally smooth function superimposed on a set of piecewise 
constant functions. Extensions of our method to more general piecewise smooth 
functions will be reported at a later date, as will its application to magnetic 
resonance imaging. 



A Proof of Proposition 13.11 

To prove Proposition 13.11 wc will require some results concerning the asymp- 
totics of integrals of the form (|3.1|) near a point of stationary phase. 



Lemma A.l Consider a curve "/{t), proceeding at unit speed, for t G [0,L]. 
Suppose k ■ 7'(0) = 0, A; • j'{t) ^ /or t € (0, L\, and k ■ j"{t) does not vanish. 
Let kg = k/ k 
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Then 



1 



i \k\ Jo 



1 



1^1 



3/2 



„ifc-7(0) 



. ' -7(1/2, -I \k\ + ^ 

2(A:e • 7"(0)/2)i/2 ' ' '^^ |fc|2 



1 



„ifc-7(0). 



-1/2 



2{ke ■ 7"(0)/2)i/2 

i?(fc) 



r(i/2) 



|fc|-3/2,....(o)^-^^:r(i/2,-.|fc|/3) 

2(fcfl-7"(0)/2)i/2 



T/ie remainder R{k) is hounded by: 

Rik)\<l + l 



ke ■ i"{t) 




k^ ■ 7'(i) 1 




ke-YiO 


fee • 7'(i) 


fce-y(i) fcf,-7"(0)i 



(A.lb) 



Here, 7(0, z) and r(a, z) are </ie incomplete Gamma Junctions jl, Chapt. 6.5, 
p.p. 269] and (3^d- j{L) - d ■ 7(0). 

Proof. This is a standard application of stationary phase [19l Section 3.13]. 
Recahing that , kg = k/ \k\, the unit direction of k, let us define the variable 

V = kg-"f[t) — k0-"f{Q) and (i = kg-j{L) — kg--f{0), so that dv = k0-"f'{t)dt. Since 

V = kg ■ 7(t) — kg ■ 7(0) = kg ■ 7"(0)i^/2 + 0{t^), a straightforward calculation 
shows that 



fce-7'(iM) 2(fce-7"(0)/2)i/2 



as w — > 0. 



Thus, 



7o 7o /fee •7'(iW) 



\k\ 



■i \k\v 



2(A;e-7"(0)/2)i/2 



-dv 



\k\ 



-1/2 



|fc| 



\kg-j'itiv)) 2(fce-7"(0)/2!)i/2^ 

ji/2 



1/2 



|fc| 



2(fc9 • 7"(0)/2)i/2 
^Aklv ( ki ■ i{t{v)) 



7(1/2, -*|fc|/3) 



-1/2 



A:e-7'(t(")) 2(fce • 7"(0)/2!)i/2 



dv (A.2) 
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We must now bound the remainder, the last hne of (jA.2[) . This can be done via 
Theorem 12.3 of [19l p.p. 99], which states that the integral is bounded by the 
total variation norm of the integrand plus the value at the endpoints, i.e.: 



-ze 



dv 



with 



le 



.fe..(o)^ fc^7"(o) 



ke ■ 7"(0) 



\E\ < 



1 



ke-Y{t) kg-Y'{0)t^ 



TV 



To compute the total variation norm, first use the expansion: 



E (A.3) 



d-Y{t) ke-Y'{0)t 



kg ■ Yit) 



1 



kg-j"{0)t 
1 



i + (Vi -(fcg-yW)^-i) 

kg ■ Y{t) 



kg ■ i'{Q)t - kg ■ Y{t) ^l-{ke- l'{tW - 1 



kg ■ 7"(0)t 



kg-Y{t)kg-Y'{0)t 



kg ■ Y{t) 



(A.4) 



[0,/?]. 

It is an exercise in elementary calculus to show that ||/(^)|| 



([0,1]) 



1 where 



f{z) — (Vl — — l)/z\ this, combined with the fact that kg ■ Y{t) G [0, 1] and 
kg •7'(i) is monotonically increasing, shows that the last term of (|A.4p has total 
variation less than 1. 

To bound the first term, we begin by using Taylor's theorem: 

kg ■ i(t) - kg ■ Y'(0)t = R2{t) = / {l/2)kg ■ Y" {z){z - t)dz 



Define R{t) = t ^R2{t). Then we can write: 



kg ■ Y'{0)t - kg ■ Y{t) 



kg-Y{t)kg-Y'{0)t 



< 



TV 



d kg-Y'{Q)t-ke-Y{t) 



dt kg-Y{t)ke-Y'W 



1 



kg ■ 7"(0) 



d R{t) 

dtkg-Yit)' 



-dt = 



1 



kg-Y'iO) 



kg-Y{t)R'{t)-R{t)kg-Y'{t) 



{ke ■ Y{t)Y 



dt 



(A.5) 
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Note that: ^ 

R'{t) = - ( (l/2)fce • -i"'{z)z/t^dz 
Jo 

Substituting (|A.6|) into (|A.5P yields: 



(A.6) 



ke ■ 7"(0) 



fee -7 [z) 



1 



L ft 



ke ■ 7"(0) 
1 



'0 JQ 



< 



ke ■ 7"(0) 
ke ■ i"{t) 



ke-l"'{z\ 
ke ■ l"'{z) 



{ke-l'{tW 
-ke-i{t)z-{zt-e)ke-i'{t) 



dzdt 



Jz 



{ke ■ i{t))H^ 
-ke-i{t)z-{zt-e)ke-i'{t) 



< 



ke ■ 7"(0) 
ke-Y"{t) 



ci: 



(ke ■ Y{t))H^ 
z zke-l"{t) ^ ke-i'{t) 



ke ■ Y{t)t^ {ke ■ -i'[t))H [ke ■ i'{t)f 



dzdt 



dtdz 



dtdz 



ke-l"{0) 



ci: 



z \ , ke-i'{t) 



< 



ke-i"{t) 



ke ■ 7"(0) 



z 1 



dtdz 



ke-^'{L)L ke-Y{z)z ke-Y{L) kg-Y{z) 



dz 



ke ■ 7"'(t) 



< 



ke ■ 7"(0) 
ke ■ i"{t) 



ke ■ 7"(0) 



ke-YiL)L ke-YiL) 



L L 

+ 



-dz 



2ke-^'{L) ke-Y{L)^ 



3 

< - 

- 4 



k 


■ i"{t) 






^e-7"(0 


ke ■ 1'{L) 



Thus, we have the bound that: 

^fce-7'W ke-rW) 
Tliis impUes that: 



< 1 



TV 



ke-i"{t) 






ke ■ 7"(0 


ke ■ i{L) 



(TOi = -i 



Ak 



,(L) / ki-i{L) 



ke-Y{L) ke-Y'{0)L, 



V ke ■ 7"(0) 



(A.7) 

(A.8) 

+ E 
(A.9) 
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with 



1^1 < 1 



ke ■ i"{t) 






kg ■ 7"(0 


ke-i{L) 



Since d is the normal vector to ^{t) at i = 0, and k-^ is tangent to it, we find 
that kj- ■ 7"(0) = (using the fact that 7"(0) = k(0)A^(0)). Thus, we find that: 



e'fc-7Wfci ^ |fc|i/2 



,1/2 



- le 



with 



ifc-7(0) ( 21. 



fc^7"(o) 

fce-7"(0). 



+ ^;' = |fc| 



2(fce-7"(0)/2)i/2 

,1/2 



1/2 



2(A:e-7"(0)/2)i/2 



7(1/2, -z|fc|/3) 

7(1/2, |fc| /3) + 



l^'l < 1 + T 





« • i"{t) 




ki ■ i{L) 1 




ke ■ 7"(0 


ke ■ l'{L) ^ 


ke ■ i'{L) ke ■ 7"(0)L 



Multiplying (jA.lOp by l/(i |fc|^) yields the result we seek. 
We also need the following geometric result: 



(A.IO) 



(A.ll) 



□ 



Lemma A. 2 Consider two normal vectors fci and k2 on a curve jj(t) with an 
angle 9 between them. Then: 



e/K> \tj(ki) ~ tj(k2)\ > \9\/R 



(A.12) 



Proof. The angle changes most quickly if 7^ [t) is a circle with minimal radius 
of curvature. The arc length along such a curve is \0\R > \d\/R. The arc 
length changes most quickly (w.r.t angle) along a circle with maximal radius of 
curvature, i.e. a circle of radius 1/k. □ 

The proof of Proposition 13.11 basically requires us to applying Lemma lA.ll 
to (IXTl) . 

Proof of Proposition 13.11 Note that k ■ ^'j{tj{k)) — 0. Thus, we may 
apply Lemma lA.ll to the curve ^{t) — Jj{tj{k) + t) along the segment t G 
[0, ifc(fc^) — tj{k)] and similarly to j{t) — "(j{tj{k) — t). This will give us an 
expansion over the sections of the curve where k ■ Nj{t) > 0. Repeating the 
analysis centered at tj{—k) yields an expansion over sections of the curve where 
k ■ Nj{t) < 0. 
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Applying Lemma [A . 1 1 directly yields i 

' 'fc-^'^fc^.ymdt ^ e^-^^o)^ — ^ r(i/2) 

_ 3/2^^k.^(o) ^-'^'r{l/2-^\k\(3,) 
2{ke-l"{h{k))/2Y/^ 



i\k^ Jk-Nj(t)>o 



the following: 

,•-1/2 
.gjfc-7(0). * 



j_ ifc-7(0). 



! r(i/2) 

\k\'^'' 2{ke--r{tj{k))l2Yn 
R2{k) _ 3/2 ^.fc.7(0) ^-^/^r(l/2,-z|fc|/32) 
|fc|' 2(^e--7"(ij 



= g»fc-7(tj(fe)) 



2(^,._y'(t^.(fc))/2)i/2 



|fe|-3/2^.fc.7(0) »"^/'r(l/2,-»|fc|/3i) 

2(fc0.7"(i,(fc))/2)i/2 
Here, the remainders Ri.2{k) are bounded by: 
i?i(fc) 



Ri[k) + R2{k) 
\kt 



\kf'^K,{t,{k)) \k\ 

»-V2r(l/2,-»|fc|/3i) _^ 3/2 ^.fc.7fo) »'^^'r(l/2,-z|fc|/32) 
2(fc0.7"(i,(fc))/2)i/2 2(^e-7"(ijW)/2)i/2 

(A.] 



<1 + ? 









ke-l"{t,{k)) 


fc,.7'(i,(fc^)) 



kj-iUk^)) 



I V" // _ 

ke ■ Y{t,{k^)) ke ■ r{t,{k)mk^) - t,{k)) 



(A. 14a) 



Ri{k) 



< 1 



ke-l';'{t) 




. {t,{k)-t,{-k^)) 


ko-l"{t,{k)) 


ke-i'{t,{-k^)) 



ke ■ i'{t,{~k^)) ke ■ Y'{t,{k)){t,{k) - t,{-k^)) 



(A.14b) 



Since k ■ ^'j{tj(k)) = and 7j(t) proceeds with unit speed, we find that ke ■ 
-y'Jitjik)) = K(tj{k)). Note also that kj- ■ i^{tj{±k^)) = since i^{tj{±k^)) \\ 

±k^ (and similarly ke ■ ^' {tj{—k^)) = 1. Substituting this into (jA.14p . as well 
as bounding the curvature below by k and adding them up yields: 



Ri{k) + i?2(fc) 



< 2 



3ii7rwiLo.t,(fc^)-t,(-fc^) 



^,{tjik)){t,{k^)-t,{k)) 



K,{tj{k)){tj{k) - tj{-k^)) 



(A.15) 
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Applying Lemma [A . 2 1 shows that {tj{k^) — tj{k)) > 7r/2K and similarly {tj{k) — 
tj{—k^)) > 7t/2R. Substituting this into (|A.15p and observing that Kj{tj{k)) > 
K yields: 



Riik) + R2{k) 
We must also bound 



< 2 



2 K 



Ak 

— itjik^) - tji-k^)) + — (A.16) 

TTK 



1^1-3/2 gifc-7(0) * -3/2 gifc.7(0) * 

2{ke-r{tj{k))/2y/^ 2{ke-i"{t,{k))/2Yn 



To do this, note that f3i — kg ■ ^j{tj{k^)) — kg ■ ^j{tj{k)) < —1/k (and similarly 
02 < —1/k. Thus, — |fc|/3i,2 > \kr\ K- This can be seen easily by considering 
circle tangent to Jj{tj{k)) of radius 1/R. Using Eq. (6.5.32) from [I] p.p. 263], 
combined with the estimate on the remainder stated immediately after Eq. 
(6.5.32), we observe that |r(l/2, |fc| /3jj)| < Ifc/Kp^/^. Thus, we find that: 



\k\ 



^3/2^a..7(o)£^^!r(iA-H^ ^ l^|-3/2^.fe.,(o)£^^!r(iA-H^ 



2(fc..7"(t,(fc))/2)V2 



2(fc,.7"(t,(fc))/2)i/2 



< 



\k\- 



3/2 jk-n)-'/' ^ .,.-3/2 



2{k/2) 



\k\- 



2{k/2) 



1 /2S 
< — A.17 

-{kn^ 



Repeating this analysis for the part of the curve where k ■ Nj (t) < yields the 
following: 



hAk) 



where 



Ej{k) 



\k\ 



3/2 



(A.18) 



Ej{k) 



<4 + 3- 



|7'"W| 



-^arclength(7j) 



8k 



(A.19) 



Adding this result up over j — Q . . . M — 1 and bounding pj by p yields the 
desired result. □ 



B Proof of Theorem [4^3]: Leading Order Asymp- 
totics 

We first show that the directional filters behave properly when applied to the 
leading order asymptotic terms of p{k) . 
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For this, we need to prove two facts: a) that directional filters, after being 
applied to the image, decay away from the points "fjitj{9)), with 9 the direction 
of the filter and b) that the filters yield spikes at or near the points Jj{tj{0)). 

The basis for our calculation is the following Lemma, which allows us to write 
the directional filter applied to the leading order term of (|3.2p as an integral 
over the curve 7j (t) . 

We consider only a directional filter oriented in the direction 9 ~ Q. Results 
for other directions can be obtained by rotation. 

Lemma B.l LetV{kg) be supported on the interval [—a, a], letW{kr) G L^{W, dk^) 
and let W{kr)/kr € L'^(R,dkr). Then: 

e^'''"' — , ^ V{kg)kr^'^W{kr)dk0krdkr 



\kf'^ 



Kj{tj{k)) 



a /"OO 



a Jo 



y/Kj{tj{k0)) 

^ ^^{ke)•W{N,(t,(ke)) ■ Yl,it,(ke)) - x])dke 



a y/Kj{tj{kg)) 



IT / ' V{ke{t))J^)mN,{t) ■ [7,(0 - x])dt (B.l) 



Proof. 



^~ — , ^ V{ke)kr^^^W{kr)dkgkrdkr 



\k[ 



a poo 



a Jo 



Note that the inner integral of the last line of (|B.2[) is merely the inverse Fourier 
transform of >V(fcr) evaluated at the point r = Nj{tj{kg))-[-fj{tj{k0)) — x]. Thus, 

(ESI = / , V{ke)yV{N,{t,{kg)) ■ [j,{t,{kg)) - x])dkg 

J-a K jit J [kg)) 

- / ^ V{kg{t))J^)W{N,{t) ■ [7,(i) - x])dt, 

Jt,{-a) ^ 

completing the proof. □ 



Proposition B.2 Let V{kg) be supported on the interval [—a, a], let W(fcr) € 
L'^{R,dkr) and also W{kr)IK e L^{M.,dkr). 
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Define the smallest normal and tangent distances (D^ and Dt respectively) 



Dn = , inf ^N,it)-[-f,it)-x] 

te[t,(-a),t,(a)] ■' 



(B.3a) 
(B.3b) 



Then we have the following bound on the action of the filter in the normal 
directions: 



-ik-x 



J -a 



\k\ 



3/2 



^=V{kg)kr^'^W{kr)dkgkrdkr 



Kj{tj{k)) 



< 



-||V(fce)ILi(si,,fe,)W(D^) (B.4a) 



We also have a weaker bound in the tangential direction: 



^^''^•^ Y V{k0)kr^^^W{kr)dkgkrdk 



J -a 



\k\ 



3/2 



nj{tj{k)) 



< 



2V^\\W{kr)/krhi^M,dk^^ f\\V'ike)hl l|V(fc9)|lil ||7-"W| 



,1/2 



V^\\V{ke)hi \\yV{z){z + l)\ 



Finally, 



„2 r 



Dn > cos{a) |7V,(t,(0)) • (7(^,(0)) - ^)\ - — 



2 k3 



6 



(B.4b) 



(B.5) 



Proof. 

The result (jB.4aP follows from the second to last line of (jB.l|l and the fact 
that 

\WiN,it,ikg)) ■ [j,it,ike)) - x])\ < W{Dm). 

To prove (jB.Sp . we must bound Nj{tj{kg)) ■ ['y{tj{kg)) — x]. Let to = ^j(O), 
and consider the Taylor expansion (to second order) of jj {t) . 



7" (to) 

Ijito) + l'j{to){t ~ to) + - to? + remainder 



By Taylor's theorem, the remainder is bounded by \k ■ R{t)\ < |l7"'(t)|l^o 
The first order term is: 



6Kj{t) ■ 



\Njikg) ■ 7'(to)| = |cos(fce + 7r/2)| < \sm{ke - 0)| < Ifc^l < 



35 



We use also the fact that |i — to| < a/ k. Thus, we obtain: 

\N,{t) ■ (7W -x)- N,{t) ■ (7(to) -x)\<- + --+ \t - tor 

Hi Zi Hi D 

2^1 , h"'mL^a^ 



Note also that Nj{t) ■ (7(^0) — a;) > cos{a)Nj{to) ■ (7(^0) — x). We therefore find 
that: 



\N,{t) ■ [^{t) - x]\ > 



cos(a)|iV,(io)-(7(io)-a;)| 



a 



K / 6 



(B.7) 



Taking the inf over the angles in [—a, a] yields the result we seek. 

For (|R4b)) . note that dtNj{t) ■ [-fj{t) ~ x] ^ Kj{t)j'^{t) ■ [yj{t) - x]. We can 



then multiply and divide the integrand of (|B.ip by this and then integrate by 
parts to obtain: 



tj (-a) V «i (t) ■ n w - ^\ 



V(fceW)«;;-(t) , V(fce(i))(7"W-[7,W-^]+7'(i)-7'W)^ 



2(«, W)^/^7' W • [7, W - ^] W(7'(i) • [7, W - 



\dt 



(-a) V V^j(i)7j (^) • [ij it) - x] 

V{ke{t)y^(t) V{ke{t)) {n,{t)N,{t) ■ [^,{t) -x] + l) 



2{^,{t)mAt) ■ [7, W - x] n,{t){iAt) ■ [7, W - x] 



yt 

(Bi 
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We can further simplify this to: 



V{kg{t)y^{t) V{ke{t)){N,{t)-[-t,{t)-x\ + l) 



\dt 



2{n,it)r/W,{t) ■ [lj{t) - x] ■ [7,(t) - x]y 

W{Nj{tj{ke)) ■ [l{tj{ke)) - x]) x — - 
V{ke{t,{ke)))nr{t,{ke)) 



2{K,{t,{keW/^i^{t,{kg)) ■ h,{t,{ke)) - x] 
V{ke{t,{ke))) {N,{tj{ke)) ■ bj{tj{ke)) -x\ + l)\ dke 



(B.9) 



Note that F{z) = J" f{z')dz' = / e~''^-^W{kr){ikry^dkr, and thus \F{z)\ < 
\\yV{kr)/kr\\^if^^iy Noting also that -/j"{t) ■ Nj{t) = (differentiate the 

formula 7"(t) = Kj{t)Nj{t), use the Frenet-Serret formula and dot product with 
Nj{t)), we find that < |7"'(t)|. Combining these facts and using the 

definition of Dt, we obtain: 

\mM\ 

+ F{N,{tAke))ii{t,{ke))~x]){N,Uke))i^{t^^^^^ 

( V{ke{t,{ke))) ^V{ke{t,{ke)))\\lT{t)\\^A 



< / ||>V(A:,)/fc,,|| 



< 



+ \\F{z){z + 1)||^^ jj:^ dkg 

2V^||W(fcr)/fc.|L.(R,,fc^) \\V'{kg)\\ 



T 



, x/^|lV(fce)|L, ||F(z)(^ + l)|L 



(B.IO) 



This yields (|B.4bp . and completes the proof of Proposition IB. 21 □ 

Lemma B.3 Let"/{t) be a curve moving at unit speed and having non-vanishing 
curvature, with unit tangent T{t) and normal N(t). Then: 

|7V(t).(7W-7(io))l<^ (B.ll) 
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Here, e{t) = T{t)ZT{to). 

Proof. Letting d{t) be the angle of the tangent, we find that 9'{t) = K{t). Note 
that non-vanishing curvature imphes that t{9) and 9(t) are both functions (at 
least for small t and 6). Note first that: 

Integrating this with respect to 6 shows that |7(t) — 7(to) I < Now compute: 

= -l^i-nit)T{t)) ■ (7(t) - 7(io)) + -^,N{t) ■ T{t) ^ T{t) ■ {^{t) - 7(^0)) 

This implies that: 
d 



^^iv(i)- (7(0-7(^0)) 



<|r(t).(7(t)-7(M)l<- 

Integrating with respect to 6 yields the result we seek. □ 



Proposition B.4 Suppose that V(kg) is smooth and compactly supported on 
[—a, a]. Let WP{r) be a function symmetric about r = 0, and strictly positive on 
the interval [— a^/2/«, a^/2K]. Suppose also thatW{kr) satisfies (|4.1ap . Assume 
also that ()4.4ap is satisfied. 

Then at the point 7j(ij(0)), we have that: 

372 V{ke)kr^/^W{kr)dkekrdkr 

\k\ ^n,{t,{k)) 



>^/^l|V(fc«)||^.(,,^) inf W^'(r) (B.12) 



TT 

Proof. Note that: 



k ■ ilAtjik)) -x)^ krN,(t,{ke)) ■ [lj{t,{ke)) - x] (B.13) 
Following the calculations of (jB.ip and using (|4.1bp as well as (|B.13I) we find: 

^ 372 V{kg)kr^'^W{kr)dkgkrdkr 

\k\ ^n,{t,{k)) 

-.V{ke)W{N,{tj{kg)) ■ hjit,{kg)) - X])dkg 



' —^==Vikg)g(N,{t,{ke)) ■ hj{t,{ke)) 

. ^ _ ,|)*, ,B,4) 



38 



For X — Jj(tj{0)) and a G [—a, a], Lemma lB.31 implies that: 

\N,it)-i^it)-x)\<^_<^ 



2k - 4 



where the last line follows from Lemma [i3.31 This implies that cos{Nj (t) ■ {j{t) — x 
1/a/2 and therefore: 



(jHH > / —^==V{ke)m^dke 



This is what we wanted to prove. □ 



B.l Proof of Theorem 14.31 Putting it together 

We now consider the behavior of the filter applied to the entire asymptotic 
expansion of the image. 

Lemma B.5 Suppose that yV'(fcr) satisfies (j4.4b|l . In that case 



e,a- 



<e:(W,V,a) inf (B.15) 

te[t,(-a),tj(a)] \^j{t) - X\ 



^K,{t,{k)) 

where €{W,V,a) is given by (j4.4d|) . 

Proof. As before, we study the case when 9 = 0, since the rest can be treated 
by rotation. 

Recalling (|B.5p . we define the quantities: 

D+ = sup inf \N,{s) ■ [j,{t) - x]\ 

se [tj i~a),tj (a)] t6 [*j (-")■* j (")] 

D+ = sup inf ms)-[j,{t)~x]\. 

Obviously > Djq and > Dt- Let N = Nj{s) be the normal at which 
is achieved, and T = Tj{s) be the tangent at which is achieved. Let 
p = ~ X. Since the angle between N and T is at most 2a, we find that 
l^l^ < sec(2a)^ \x\\n t)' where | • j,^ is the norm taken in the coordinate 
system {N,T). Thus, we find: 

\p\^ < sec(2a)2 j,) < 2 sec(2a)2 max{|7V ■p\,\T-p\y 
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Thus: 



sup \jj{t)-x\ 

tG[tj(-a),tj(a)] 

< y2sec(2a) sup ma-x{\N ■ {jj{t) - x)\ ,\T ■ {-fj{t) - x)\y 
= \/2sec(2Q!)max{ sup \N ■ {jj{t) ~ x)\ , sup \T ■ {jj{t) - x)\} 

= \/2sec(2a) max{i:)+, D+j 



This impUes that: 
inf 



1 cos(2a) . r 1 1 ] cos(2a) . f 1 1 
; — r > ^ — min < — r , — r f" > ^ — min < — — , -— 

l7,(i)-^|- V2 [DyD+j- V2 [Dr^'DT 



Now, substituting (j4.4bp into (|B.4a|) we find: 



,3/2 



Using (|B.4b|) . we find: 



< 



e:wv/V^llv(fce)|lii(si^dfc«) 



D 



N 



(B.16) 



(B.17) 



,|3/2 



< 



2V^||>V(fcO/fcr|Li(R,dfc^) /||V'(fce)||^i l|V(fce)|lii ||7-"W|Lo 



, 0F||V(fcg)IL. ||W(z)(z + l)| 



2«;5/2 



< 



2||W(fc,,)/fc.||^,(R,,,^ 



|U"'(t)|| 



2^2 



+ Vk||W(z)(z + 1)| 



(B.18) 
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Taking the min of (|RT7)) and (|RT8)1 . and using (|RT6)1 wc find: 



aik-yj{tj{k)) 



\k\ 



3/2 



. I 1 

< mm < -r- — , -— 



X maxjCw II V(A:e)|lLi(si,dfce) 



£1 -T II v^wyiiLi ^ — 



1v'(MIIli + I|v(MILi 



+ ^||W(z)(z + l)||^„ jj 
<(inf. ' ^ ^ 



\jj{t) — x\J y/KCOs(2a) 



max|G:w||V(/ce)|lLi(si,dfe,), 



2||W(A;,)/fc.||^,(R,,_, 



v'(fc«)ILi + ||v(MII 



2^2 



+ V^||W(2)(z + l)| 



(B.19) 



This is what we wanted to show. 

Recalhng Definition 14.21 we have the following result: 



□ 



Proposition B.6 Let x be a point not on any of the curves ^i{t). Then (with 
£(yV,V,Q;) as in Lemma \B.5\) : 



De.a 



3=0 



< 



l|v(fc«)|| 



(B.20) 



where ^geo is defined in \3. 3I\) . 
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Proof. Begin by using (|3.2p : 



A/-1 
J=0 



< 



M-1 
3=0 



^ik-,{tj{k)) 



De.a 



m 



M-l 

J=0 



De.a- 



^ik--){tj{k)) 



\k\ 



3/2 



E{k) 



De,a- 



Jk-i(t,(-k)) 



3=0 



inf 



\k\' 
1 



inf 



te[t3(-a+6i),tj(Q+e)] |7j(t) - a;| te[t-i{-a+e+n),t,(a+e+7r}] {jjit) - x\_ 

E{k) 



e:(W,V,a)2Af 
< — 77 — — 7 — maxp,- 



€(W,V,a)2M_ 
< p 



To control the last term, recall the bound Cgco in p.3p on E{k): 



Do 



E{k) 



ifcr 



E{k) 



-'^■''V{ke)kr^''^W{kr)^^dk0krdkr 



< 



W{kr] 



< 



Wikr) 



Vike)E{k) 



dkgdkr 



1/2 



dkr j Cgco |V(fce)| dfcfl 



< 



l|V(fce)|| 



□ 
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Proposition B.7 We have the following estimate at the point x — lj{tj{9)). 



■N,(t)]De^c. 



A/-1 
J=0 



e:(>V,V,a)(2Af- 1) 



Li(R,d/c^) 



l">^(fc9)llLi(Si,'ifefl) 



(B.21) 



Proof. By Proposition IB. 4[ we have that: 



gifc-7j(tj(fc)) ^ 
^S,aP , , ,3/2 ; ^ ■ Ifco g [-a,a] (^9 ) 



|fc| 



nj{tj{k)) 



> 



Pl|V(fc.)IL.(,fe,) . inf W^r) 



Applying the same reasoning as in the proof of Proposition IB. 61 we find that: 



De.a 



< 



£(W,V,a)(2Af-2) 



The second triangle inequality yields the result we seek. □ 

Proof of Theorem 14.31 Proposition IB. 61 proves (|4.4cp . 

Note that Assumption [2] implies that rf(7j(<j(^^)), 2le,a"') > S. Substituting 
this into Proposition IB . 71 proves ()4.4ep . and thus completes the proof of Theorem 
[431 □ 
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